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We introduce a stochastic algorithm that acts as a prime number generator. The dynamics of such 
algorithm gives rise to a continuous phase transition which separates a phase where the algorithm 
is able to reduce a whole set of integers into primes and a phase where the system reaches a frozen 
state with low prime density. We present both numerical simulations and an analytical approach in 
terms of an annealed approximation, by means of which the data are collapsed. A critical slowing 
down phenomenon is also outlined. 
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From the celebrated coincidence in 1972 between H. 
Montgomery's work on the statistics of the spacings be- 
tween zeta zeros and F. Dyson's analogous work on eigen- 
values of random matrices, we have seen, somewhat un- 
expectedly, how number theory and physics have built 
bridges between each other. These connections range 
from the reinterpretation of the Riemann zeta function as 
a partition function [l| or the focus of the Riemann Hy- 
pothesis via quantum chaos 0, to multifractality in the 
distribution of primes @ or computationalphase transi- 
tions in the number partitioning problem to cite but 
a few (see Js'] for an extensive bibhography) . 
Prime numbers are mostly found using the classical sieve 
of Eratosthenes and its recent improvements [6]. Ad- 
ditionally, several methods able to generate probable 
primes have been put forward In this paper we 

study a somewhat different algorithm from those men- 
tioned above, based on an artificial chemistry model in- 
troduced by Dittrich Q that generates primes by means 
of a stochastic integer decomposition. 
Suppose a pool of positive integers {2, 3, Af}, from 
which a subset of numbers is randomly extracted. No- 
tice that the number 1 is ignored and that repetitions are 
allowed in the subset. Given two integers a and h (taken 
from the subset of A^ numbers), the reaction rules of the 
algorithm are: 

• Rule 1: li a — b then no reaction takes place, and 
the numbers are not modified. 

• Rule 2: If the numbers are different, say a > 6, a 
reaction will take place only if 6 is a divisor of a, 
i.e. if there exists an integer c — a/b. Then a is 
eliminated from the subset and substituted by c. 

• Rule 3: On the other hand, if a is not divisible by 
6, then no reaction takes place. 

The stochastic algorithm goes as follows: after choos- 
ing a subset of A^ numbers from the pool {2,3,...,Af} 
two numbers a, h belonging to that subset are picked at 
random; then the reaction rules are applied to them. We 
consider N repetitions of this process as one time step in 



order to have a parallel updating. Notice that the posi- 
tive reactions tend to decompose numbers, thereby this 
process when iterated may generate primes. We say that 
the system has reached the steady state when no more 
reactions are achieved, either because every number has 
become a prime or because rule 2 cannot be satisfied any- 
more: the algorithm then stops. 

In what follows we will firstly present the phase transi- 
tion that the system seems to exhibit. Second, we will try 
to interpret this phase transition in terms of a dynamical 
process embedded in a directed catalytic network, intro- 
ducing subsequently a proper order parameter. Some 
analytical arguments in terms of an annealed approxi- 
mation will then be outlined in a third part, where a 
data collapse is provided. Finally, a critical slowing down 
(easy-hard-easy pattern in the algorithmic phase transi- 
tions language) is pointed out. 

A preliminary indicator of the system's behavior may 
be the unit percentage or ratio of primes r, that the sys- 
tem reaches in the steady state. This parameter will char- 
acterize the ability of the algorithm to produce primes, 
for a given TV. In figure [1] we plot the behavior of r 
versus A^ for several pool sizes M . We clearly see that 
two separated regimes arise: the first one is characterized 
by small ratios (low proportion of primes in the station- 
ary state) while in the second one every single number 
of the system will end up as a prime. The system thus 
exhibits a sort of phase transition. Note that r is not a 
well defined order parameter since it does not vanish in 
the disordered phase. This is due to the fact that fol- 
lowing the prime number theorem Q, in a pool of size 
M there are typically M/ ln(M) primes. This residual 
value of r is not related to the algorithm dynamics. In 
fact, when A^ is small the number of reactions until the 
system reaches the steady state is quite small. Therefore, 
the residual ratio r « 1/ ln(M) is the relevant contribu- 
tion p^ . It comes thus necessary to define an adequate 
order parameter that would properly describe the former 
phase transition. 

Let us now see how this phase transition can be un- 
derstood as a dynamical process embedded in a catalytic 
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FIG. 1: The ratio r of prime numbers in the steady state 
as a function of A'^ for different pool sizes M. Each point is 
the average of 2 x 10* reahzations. For small values of A^, 
the value of r is only related to the amount expected from a 
random sample. For large values of A'^, r tends to one: the 
algorithm is able to reduce the whole system into primes. 
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tribution is P{k) ^ with A = 2. In our system, fixing 
N is equivalent to selecting a random subset of nodes in 
this network. If a and b are selected they may react giv- 
ing a/b = c; in terms of the network this means that 
the path between nodes a and b is traveled thanks to the 
catalytic presence of c. We may say that our network is 
indeed a catalytic one [l^, [l^ where there are no cycles 
as attractors but two different stationary phases: (i) for 
large values of TV all resulting paths sink into prime num- 
bers, and (ii) if N is small only a few paths are traveled 
and no primes are reached. Notice that in this network 
representation, primes are the only nodes that have input 
links but no output links (by definition, a prime number 
is only divisible by the unit and by itself, acting as an 
absorbing node of the dynamics). When the temporal 
evolution of this algorithm is explored for small values of 
N , we observe that the steady state is reached very fast. 
As a consequence, there are only a few traveled paths 
over the network and since N is small the probability of 
catalysis is small as well, hence the paths ending in prime 
nodes are not traveled. We say in this case that the sys- 
tem freezes in a disordered state. In contrast when N is 
large enough, many reactions take place and the network 
is traveled at large. Under these circumstances, an arbi- 
trary node may be catalyzed by a large N — \ quantity 
of numbers, its probability of reaction being high. Thus, 
in average all numbers can follow network paths towards 
the prime nodes: we say that the system reaches an or- 
dered state. 

In the light of the preceding arguments, it is meaningful 
to define an order parameter as the probability P{N, M) 
that the N numbers extracted from M be primes once 
the stationary state is reached. In figure [5] the relation 
between the order parameter P and the control parame- 
ter N related to the same simulations that in figure [T] is 
depicted. Note that P is now a well defined order param- 
eter, as opposed to r. In each case, Nc{M) is the critical 
value separating the phases P = and P 0. Observe 
in figure m that Nc increases with the pool size M. In 
order to describe this size dependence, we need to find 
some analytical argument by means of which to define a 
system's characteristic size. As we will see below, this 
one will not be M as one would expect. 



FIG. 2: Order parameter P, defined as the probability that 
all numbers are prime in the steady state versus the control 
parameter A^, for the same pool sizes as for figure [1] (simula- 
tions are averaged over 2 x 10* realizations). Inset: scaling of 
versus M/ln(M), for pool sizes M = 2^° ,2^^ , ...,2^\ 



network having integer numbers as the nodes. Consider 
two numbers of that network, say a and b {a > b). These 
numbers are connected (a b) if they are exactly di- 
visible, that is to say, if a/6 = c with c being an inte- 
ger. The topology of similar networks has been studied 
in [ll|, [l2, concretely in [l^ it is shown that this 
network exhibits scale-free topology the degree dis- 



Note that non-trivial correlations between the values of 
the N numbers take place at each algorithm's time step. 
This leads to highly complex, analytically untractable 
dynamics. We can try however an annealed approxima- 
tion in order to break these correlations, assuming that 
at each time step, the N numbers are randomly gener- 
ated. In figure [3] we depict for different values of M , 
the resulting simulations of the function 1 — where 
q — q{N, M) is the probability that no pair of randomly 
chosen N numbers from M be divisible between them. 
Notice that once M is fixed, there is a certain value of 
N from which the probability of finding at least one re- 
acting pair is almost 1. The behavior of 1 — g follows 
qualitatively the behavior of the order parameter P. In 
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percolation processes, since the choice of the percolation 
threshold, related to the definition of a spanning cluster, 
is somewhat arbitrary in finite size systems [Ol- After 
some algebra and taking a leading-order approximation, 
we find the scaling relationship: 
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In the inset of figure [3] we plot in log- log the scaling be- 
tween Nc and M/ ln(M) in the annealed system, which 
follows equation ([3]) with a = 0.48 ± 0.2 (note that 
within the error bar, there is indeed independence be- 
tween pairs). The goodness of the former scaling sug- 
gests that the above ansatz is acceptable. 
In appearance, in the prime number generator system. 



FIG. 3: Probability 1 — q{N, M) that at least two numbers 
from be exactly divisible between them, for different values 
of M (from left to right: 2", 2", 2l^ 2" and 2^"^). Each point 
represents the average of 10^ realizations. Inset: scaling, in 
the annealed system, of Nc (defined such that q{Nc, M) — 0.5) 
versus M/ln(M) for different values of M = 2", 2^^ 



fact, this annealed approximation suggests that once M 
is fixed in the algorithm, from a certain N we can be 
sure that at least one reaction will take place. As long as 
reactions produce new numbers while the total number 
N is conserved, reactions will then take place until the 
system reaches a stationary state of only primes. 
The probability p{M) of a reaction between two ran- 
domly chosen numbers from the pool M, that is to say, 
the probability that two numbers from {2,3,...,Af} be 
divisible is: 
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Obviously, I — p{M) is the probability that two ran- 
domly chosen numbers from {2,3, ...,M} are not divis- 
ible. From a set containing TV randomly chosen num- 
bers, the N{N — I)/2 different pairs that we can form 
are not independent, therefore the probability q{N, M) 
is not simply (I — p{M))-'^^^~^'>/^ . Correlations between 
pairs can however be taken into account through the fol- 
lowing ansatz: 
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where the exponent a characterizes the degree of de- 
pendence between pairs. For convenience, we assume 
that the threshold Nc{M) in this annealed approxima- 
tion is the one for which half of the configurations reach 
an ordered state, that is to say, the values for which 
q{Nc,M) = 0.5. This procedure is usual for instance in 
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FIG. 4: Collapse of the order parameter P for different values 
of M = 2^^, 2l^ 2l^ 2^"^, with a = 0.59. 

the characteristic size is M, however the annealed ap- 
proximation suggests that the true characterization is 
M/ ln(M). From the point of view of the network this is 
very reasonable since the amount of primes that we can 
reach increases with M in a non-linear trend, in fact it 
grows asymptotically as M/ ln(M) Q . Coming back to 
the prime generator system, in order to prove our forego- 
ing conjecture, in the inset of figure [2] we plot in log-log 
the values of Nc versus M/ln(Af). The scaling suggests 
the same relationship as equation ^ with a scaling ex- 
ponent a = 0.59 ± 0.05, that is to say, we find that the 
transition point shows critical behavior, as expected. 
In order to seek consistency, in figure|4]we collapse several 
curves P{N, M) for different pool sizes M . For that task 
we apply generic techniques of finite size scaling, where 
the size scaling is given by the function G{M / \n.{M)) — 
(M/ln(M))", with a = 0.59. Note that the data col- 
lapse is excellent. This fact gives credit to the scaling 
ansatz and provides consistency to the full development. 
As long as N is indeed an extensive variable and in order 
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FIG. 5: Characteristic time r as defined in the text ver- 
sus A'^, for different pool sizes, from left to right: M = 
2^°, 2", 2^^, 2", 2". Every simulation is averaged over 2 • 10* 
realizations. Note that for each curve, t(N, M) reaches a 
maximum in a neighborhood Nc{M). 

to find the transition point in the thermodynamic limit, 
it is meaningful to define a reduced control parameter 
n — hi(Af ) 1 which is now an intensive variable. In the 
thermodynamic limit, we find = 0. 
Some other ingredients characterizing the phase transi- 
tion can be put forward. First, we may argue that the 
responsible for the phase transition is a breaking of sym- 
metry between steady state distribution of the N inte- 
gers. As a matter of fact, we can distinguish a disordered 
phase {N << Nc) where the steady state distribution of 
the N numbers is uniform (each number appears with 



the same probability) from an ordered phase (TV >> Nc) 
where this distribution is in turn a power law . 
Second, in figure [5] we plot for different pools the be- 
havior of the characteristic time r versus N. t{N, M) is 
defined as the number of time steps per number that the 
algorithm needs to perform in order to reach the steady 
state: this parameter characterizes the relaxation time of 
the algorithm. Note that in each curve t{N, M) reaches a 
peaked maximum in a neighborhood of Nc{M) (any shift 
is due to finite size effects). Moreover, for larger pools 
this maximum is larger, diverging in the thermodynamic 
limit: this behavior is related with a critical slowing down 
phenomenon p^ . which in algorithmic phase transitions 
is known as an easy-hard-easy pattern. 

In summary, in this paper we explored a stochastic 
algorithm that works as a prime number generator. 
Many ingredients suggest the presence of a phase 
transition in the system. This unexpected behavior 
raises some interesting related questions that will be 
considered in further work, namely: how the character 
of such transition is related to with the computational 
complexity of the algorithm (l9l|? In which way the 
algorithm, which produces primes by means of stochastic 
decomposition, is related to the integer factorization 
problem and cryptography? 
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